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Abstract 

We study the dynamics of the interface between two immiscible fluids in contact with a chemically 
homogeneous moving solid plate. We consider the generic case of two fluids with any viscosity ratio 
and of a plate moving in either directions (pulled or pushed in the bath). The problem is studied 
by a combination of two models, namely an extension to finite viscosity ratio of the lubrication 
theory and a Lattice Boltzmann method. Both methods allow to resolve, in different ways, the 
viscous singularity at the triple contact between the two fluids and the wall. We find a good 
agreement between the two models particularly for small capillary numbers. When the solid plate 
moves fast enough, the entrainment of one fluid into the other one can occur. The extension of 
the lubrication model to the case of a non-zero air viscosity, as developed here, allows us to study 
the dependence of the critical capillary number for air entrainment on the other parameters in the 
problem (contact angle and viscosity ratio). 
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I. INTRODUCTION 



The problem of a solid plate pulled from a liquid bath has attracted considerable atten- 
tion in the past including the seminal contributions from Landau, Levich, Derjaguin and 
Bretherton l|-|3|. The problem has continued to be investigated, with a particular focus 
on the situation of partial wetting for which a dynamical wetting transition is observed: a 
liquid film is deposited only when a critical speed of withdrawal is exceeded (see {4, 5| for 
reviews). Much less is known about the reversed case, where the solid is plunged into the 
liquid. Again, a dynamical wetting transition has been observed, now resulting into the 
entrainment of an air film or air bubbles jsH]. Despite the viscosity contrast between the 
liquid in the reservoir and that of the surrounding air, the dynamics inside the air is very 



important for this process. The perturbation analysis by Cox [12[ suggested that the critical 
speed is inversely proportional to the viscosity of the liquid, l/r/^, with logarithmic correc- 
tions due to the viscosity of the air, rig. This is similar to the scenario for air entrainment 



by viscous cusps 



13| . such as observed for impacting liquid jets 



Q. 



Recently, in the experiment of plunging a plate into reservoirs of different liquids Marc- 
hand et al. jsl observed that the dependence on rji is much weaker than predicted; enforcing 
a power-law fit to their data would give a small exponent, in between -1/2 and -1/3 rather 
than the expected -1. This implies that air viscosity plays an important role on the on- 
set of air entrainment even if it is orders of magnitude smaller than the liquid viscosity. 
The importance of air was already highlighted in similar dip-coating experiments, where 
a reduction of the ambient pressure was shown to significantly enhance the critical speed 
of entrainment [g]. Yet, this raises another paradox: the dynamical viscosity of the air is 
virtually unchanged by a pressure reduction. 

In this paper we provide a new framework to study air entrainment by advancing contact 
lines, in which the two-phase character of the flow is taken into account. The usual lubrica- 
tion approximation is valid for small (liquid) contact angles and does not take into account 
the flow inside the gas. Since both assumptions are no longer valid near the onset of air 
entrainment, we extend the lubrication approximation such that it allows for large angles 
and a nonzero gas viscosity. This generalized lubrication theory is validated by comparing 
to Lattice Boltzmann simulations. We show that the meniscus shapes obtained from the 
generalized lubrication model agrees well with the simulations. Note, however, that such 
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Phase A 



Phase B 




FIG. 1: Schematic diagram showing a plate moving with speed Up with respect to an interface of 
two immiscible fluids. The plate has an inclination a, here drawn with a = tt/2 for the case of a 
vertical plate. The interface touches the wall at the contact line with an angle assumed to be the 
same as the equilibrium contact angle Be- The meniscus profile is described by h{x) or the local 
angle 6{s), where s is the arc length of the interface measured from the contact line. The total 
meniscus deformation is quantified by A, the distance between the contact line and the level of the 
bath. 



simulations are limited in terms of viscosity ratio and spatial resolution (i.e. the separation 
between the capillary length and the microscopic cutoff given by the interface width or, 



equivalently, by the effective slip length 



15|). and cannot achieve experimental conditions 



for air entrainment. To compare to experiments and to explore the parameter space, in 
particular the importance of air viscosity r]g, we therefore provide a detailed study using the 
generalized lubrication model. 
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II. FORMULATION 



We consider a smooth, chemically homogeneous solid plate translating across an interface 
of two immiscible fluids at a constant speed Up (positive/negative for plunging/withdrawing). 
As sketched in Fig. [H the two fluids are contained in a reservoir much larger than all the 
lengths of the problem and the plate can be inclined to any angle a. If the plate is not 
moving {Up = 0) there is no flow in the fluids, and the interface equilibrates to a static 
shape due to balance between capillarity and gravity. The interface makes an equilibrium 
angle 6^ with the solid as a result of intermolecular interaction between the three phases 
at the contact line. Neglecting the contact angle hysteresis, 6^ takes a well-defined value 
determined by Young's law. The contact line equilibrates at a distance A above the bath, 
which can be expressed as 



A = ±£^^2[1 -cos(«-^e)]- (1) 

Here, = ^Jj^^z^^ is the capillary length, defined by surface tension 7, gravity g and the 
density difference pe — pg. The ± sign depends on whether 6e is smaller (+) or larger (-) 
than the plate inclination a. When addressing the transition to air entrainment, we will 
consider the upper phase A in Fig. [T] to be gaseous, while phase B is a liquid. We therefore 
use subscripts and '(7' to indicate liquid and gas phase respectively. 

When the plate is moving, the viscous drag generated by the moving plate gradually 
deforms the fluid-fluid interface. As long as the speed of the plate is lower than a threshold 
value, the meniscus equilibrates at a new distance A from the bath level (see Fig. [1]). For 
positive Up the plate is moving downwards so that A is lower than the static equilibrium 
height, while the opposite holds for negative Up. The contact line is assumed to be straight 
so that the problem can be treated as two-dimensional. The interface is described by the 
film thickness profile h{x). The origin x = is chosen at the contact line. When the plate 
is moving beyond a critical speed, the meniscus can no longer equilibrate to a steady state. 
In the case of receding contact lines (withdrawing plate) this leads to the deposition of a 



liquid film 



16l-l22|. while air will be entrained for advancing contact lines |6-10,|23[J25|. 



In this study we only focus on viscous flows, for which the Reynolds number is assumed 
to be zero. Thus for incompressible fluids, the flow fields in the fluids are described by Stokes 
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equations and continuity, 

VgV^Ug - Vpg - = 0, V-Ug = 0, (2) 

and 

r]eV^Ui - Vpi - = 0, V ■ ue = 0, (3) 

where Ug and U£, Pg and pi, $g and are the corresponding velocity fields, pressures 
and gravitational potentials in phase A and phase B respectively. When considering air 
entrainment, phase A is assumed to be a gas of viscosity r]g, and phase B is a liquid of 
viscosity rji. The relative viscosity is expressed by the viscosity ratio R = Tjg/rji which, in 
practical situations, can be very small. 

To solve for the fiow fields and the interface shape, we need to specify appropriate bound- 
ary conditions. At the steady interface, the velocities parallel to the interface are contin- 
uous and the velocities normal to the interface vanish, so 

u\ = u\ (4) 

and 

< = K = 0. (5) 

In each phase, we define the normal stress = n - a ■ h at the interface, where n is the unit 
vector normal to the interface and the stress tensor a is defined as (in cartesian coordinates) 



, dui duj , , , 



The normal stress discontinuity across the interface is related to the curvature k and to the 
surface tension 7 by Laplace's law 

a," - < = ^K. (7) 

By contrast, the tangential stress component a* = t ■ a ■ ri is continuous across the interface, 

al = a\. (8) 

At the solid/fiuid boundary {y = 0), the velocity normal to the wall Uy{y = 0) vanishes as 
no penetration of fiuid through the solid is allowed, i.e. 

Uy{y = 0) = 0. (9) 



L Uxiy = 0)) the situation is more 



26| : imposing a no-slip boundary 



Regarding the velocity component parallel to the wa 
subtle because of the moving contact line singularity 
condition leads to diverging stress fields and calls for a microscopic mechanism of regular- 
ization. In the following section we present two methods to solve the fiow equations, which 
naturally involve different regularizations of the singularity: a generalized lubrication model 
and Lattice Boltzmann simulations. The microscopic boundary condition will therefore be 
discussed separately below. 



III. METHODS 



In this section we present two methods to determine the meniscus shape sketched in 
Fig. HJ We first present a model that can be considered as a generalization of the standard 
lubrication approximation. We then present the Lattice Boltzmann method, which is a 
rather different approach to solve for the flow and the meniscus shape. The models will be 
refereed to as GL (Generalized Lubrication model) and LB (Lattice Boltzmann). 



A. Generalized Lubrication model 



1. Derivation 



The lubrication approximation has been a very efficient framework to deal with thin film 
fiows {2^]. This systematic reduction of the Navier-Stokes equations is very suitable for nu- 
merical simulations and in many cases allows for analytical results jj]. It is usually derived 
for fiows that involve a single phase that constitutes a "thin" film, i.e. the slopes dh/dx are 
assumed small. However, the expansion parameter underlying the analysis is not the inter- 
face .ope, but tKecapm.yn„.be..Ca..W.0. ™s means that a lubrication-tj^je 
analysis can be performed whenever surface tension dominates over viscosity. Indeed, it was 



shown in 



that the lubrication approximation can be generalized to large interface angles 



6, giving a perfect agreement with the perturbation results by Voinov [29[ and Cox |12 |. 
Here we further extend this approach by taking into account, besides of the effect of a large 
slope, also of the viscous fiows on both sides of the interface. The goal is to provide a model 
that can deal with moving contact lines in cases where both phases are important (as in 
Fig. H]). In particular, this will allow us to study the air entrainment transition. Let us 
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now derive this generalized lubrication model. As mentioned in the preceding section, the 
interface curvature /t is determined by the normal stress difference across the interface. In 
curvilinear coordinates, we write 

7/^ = 7^ = < - <, (10) 

where 9 is the local interface angle and s the arc length (see Fig. [1]). The normal stresses 
have to be determined from the flows in the fluids, which themselves depend on the shape 
of the interface. For the usual lubrication theory in which the interface slope is small, the 
leading order contribution to the flow reduces to a parabolic Poiseuille flow inside the film. 
This can be generalized to two-phase flows and large interface slopes: as long as the capillary 
number is small, the interface curvature is small as well and the leading order velocity field 
is given by the flow in a wedge (Fig. |2]). The wedge solutions have been obtained analytically 



by Huh & Scriven [26|, for any viscosity ratio R = rjg/rji and for any wedge angle 6. Fig. 
[2] shows the corresponding streamlines. Based on these exact solutions, the local normal 
stress can be determined, thus giving the local curvature of the interface through Eq. fllOl) . 

We denote the quantities derived from the Huh-Scriven solutions by capital symbols, e.g. 
normal stress is denoted by S", velocity by U and pressure by P. For the Huh-Scriven 
solutions, it turns out that the non-isotropic part of the normal stress at the interface 
vanishes so that 

T/" = -P. (11) 
Approximating the normal stresses by the Huh-Scriven solutions, Eq. (ITII]) becomes 

7^ = - S," = P,- Pe. (12) 

Since Pg — Pi are defined up to an integration constant, it is convenient to differentiate Eq. 
f lT2|) once with respect to s, giving 
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d^e _ dPg dPe 
ds^ ds ds 



VPg - VP, 



int 



Cs- (13) 



The index "int" indicates that the quantities inside the brackets are to be evaluated on the 
interface and Cg = t is the unit vector tangent to the interface. 

When Stokes equation ([2]) is applied, Eq. f[T^ can be rephrased in terms of Huh-Scriven 
velocity fields 

^ V9^'% - Vi^'^Ue - V(<l', - $,)J.^^ ■ es. (14) 



^ds^ 
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The viscous contributions on the right hand side can be expressed in terms of R and 9 in 
the form 



int 



3veUpf{e,R) 
h? 



(15) 



where 



f{0,R) 



2sm'e[R^Me) + 2RMe) + fi{7r 



3[Rf,{e)f2{7r-e)-f,{7r-e)h{e)] 

f2{9) = ^ -sin^cos^ 

f,{e) ^ {9(71 -e) + sin' e). 



The gravity terms in ([Hj) can be simphfied to 



■(pe - P5)5'sin(6' - a), 



(16) 



(17) 



where pg and pi are the densities of the two phases. 

The final result of the analysis is a generalized form of the lubrication equation, which 
after scaling all lengths with the capillary length L/ = 



(fe 3Ca 



(pe-pg)g 



f{9,R) + sm{e - a). 



becomes 



(18) 



In the reminder we will use the same symbols for rescaled lengths. Once more, the capillary 
number is defined based on the viscosity of the liquid, Ca = •qiUp/'-f. This equation must be 
complemented by the geometrical relation 



dh . 
— = smO. 

ds 



(19) 



Note that for all numerical examples in the rest of the paper, we consider a vertical plate, 
for which a = n/2. 

One easily verifies that the standard lubrication equation is recovered when taking the 
limit of vanishing R, 9 and a. Namely, /(0,0) = —1 and ( IT8|) reduces to 



h"' = ^-h' + a. 



(20) 



When considering a single phase with arbitrary angle, one recovers the equation previously 



proposed in 



28 



with f{d, 0) instead of f{d, R). 



261 ] ■ and is used here to derive a 



FIG. 2: Streamlines of the flow in a wedge of angle ^«,ed, which in this case is close to vr. This 
flow solution has been derived analytically by Huh & Scriven 
generalized lubrication model. 

2. Slip boundary condition 

It is important to note that Eq. flTHl) is derived from the Huh-Scriven solution with no-slip 
boundary condition. Near the contact line, however, this induces a divergence of the pressure 
and of the shear stress, which scale as ~ rjiUp/h and ~ TjgUp/h in the liquid and the gas 
respectively. A purely hydrodynamic approach to regularize the singularity is to impose a 
slip boundary condition on the solid wall. No analytical solution for the flow in a wedge with 
slip can be obtained. However, as contact line flows are only mildly affected by the details 
of the microscopic conditions |30|, we proceed by a phenomenological regularization. We 
therefore consider the standard lubrication equation, which can be derived including a Navier 
slip boundary condition. It reads 

^ 3Ca _ ^/ .21) 



where is the slip length. In comparison to fl20|) . the effect of slip can be summarized by 
a correction factor h/{h + SA^) for the viscous term. Indeed, this weakens the singularity 



such that the equations can be integrated to /i = 



31| . We simply propose to use the same 



regularization factor for the viscous term in fllSp . i.e. 

where we have assumed the shp length to be independent of R. The appropriate boundary 
conditions are that the equihbrium contact angle 9e is recovered at the contact line, and 
that the interface attains the angle of the reservoir at infinity: 

h{s = 0) = 0; e{s = 0) = e^- e{s ^oo)=a. (23) 

The meniscus shape is now fully determined by the lubrication equation fl22|) . geometry 
f ll9p and boundary conditions ( I23p . For a given value of the capillary number Ca, the 
model parameters are the viscosity ratio R, the contact angle O^., and the microscopic length 
Xg. Below, we compute the shape of the meniscus for different parameters by numerical 
integration, using a 4th order Runge-Kutta numerical scheme. As the boundary conditions 
are imposed at different locations, the solution is determined using a shooting algorithm. 

B. Lattice Boltzmann Method 

In this section we discuss some of the key features of the LB method. In a LB model, the 
following discrete Boltzmann equation is solved for the single particle distribution function 
/i,Q,(x, t) over a 2D square lattice: 

/,,„(x + e^At, t + At)- t) = fii,„(x, t) , (24) 

where 

fii,a(x,t) = -[fi,a{^,t) - f-y:>C,t)], 



is the single time relaxation, linear BGK collision operator 32|, /f^(x, t) is the discrete 
Maxwell distribution function defined as: 



/^a(x,t) = paWi 

where 



^2 



(25) 



pa — ^ ^ fi,ai Pa^a — ^ ^ fi,a^i ; (26) 



and Wi, Bi in are the weights and the corresponding lattice speeds respectively [33|, |34 |. 



The weights, Wi, corresponding to the 2-dimensional and 9 velocity model, D2Q9, are given 
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by Wq = 4/9, ^1 = ^2 = ^3 = ^4 = 1/9, and = Wq = Wj = Wq = 1/36. The total 
fluid density is p = J2aPa and the total hydrodynamic velocity is u = J2a Pa'^a/ P- The 
effective kinematic viscosity is related to the relaxation time of the different components 
= J^aclicTaCa — 0.5) jsSi, Ca = Pa/ P is the Concentration, and Cg = l/V^ is the speed of 
sound on the lattice. In absence of an external force, each component satisfies the ideal gas 
equation of state p = c^p. For multicomponent simulation we are using two distribution 
functions {a = 1,2), whereas for multiphase simulations we restrict ourselves to only one 
distribution function (a = 1). 



1. Multiphase/multicomponent model 

The multicomponent /multiphase algorithm is based on a standard Shan-Chen lattice 
Boltzmann method 35l437l|. The non-ideal nature of the fluid is introduced by adding an 
internal force and shifting the lattice Boltzmann equilibrium velocity as 

u^q = u' + ^^, where u' = . (27) 

Pa Z^Pa/CTa 



For the non-ideal interaction the force F„ in the Shan-Chen model 



35 



37( 1 is given by: 



F„ = -GaQ'^a(x) J2 + ei)ei (28) 

where {a, a'} = {1,2} are indices for two fluid components while the coupling parameter 
Gaa' is the strength of the interaction and determines the surface tension in the model. 
This force allows for the spontaneous formation of an interface between the different fluid 
components, i.e. no interface tracking is needed. For multicomponent simulations Gu = 
G21 = G, Gn = G22 = 0, ipa = Pa- In the case of multiphase simulations a = 1, ip = ^ — 
exp(— p). The equation of state is modified to p = c^(pi + P2) + G'c^pip2 and p = c^p+ ^cj.ip^ 
respectively for multicomponent and multiphase simulations, where the first term correspond 
to the ideal gas and the second term is the non-ideal part due to the external Shan-Chen 
force. Many validation studies exist, showing that the hydrodynamical fields, u(x, t), p(x, t) 
satisfies the Navier-Stokes equations with a non-ideal Pressure tensor, under a suitable 
multiscale Chapman-Enskog expansion. 



11 



2. Boundary conditions for LB simulations 



The no -slip boundary condition for the fluid corresponds to the bounce back boundary 



condition 



33| for the distribution functions /^^^(x, t) defined at the boundary nodes. 



The surface wetting for multicomponent simulations is introduced by adding an additional 



force at the wall 



38 



^ads 



s x + e,; e,; 



(29) 



where, s(x + e^) = 1 for a wall node and is for a fluid node. The parameter G'^'^ can 
be varied to control the wetting properties of the wall; in all our simulations we have used 



ads 



G'^'^'^. Similarly for multiphase simulations we fix a wetting parameter, for the 



nodes in the wall and calculate the Shan-Chen force 



at the wall 



-Gi){p^) Wis{yi + ei)ei 



(30) 



where, s(x + e,j) = 1 for a wall node and is for a fluid node. The parameter is varied to 
control the equilibrium contact angle at the wall. Let us stress that all LB methods, inde- 
pendent of the underlying kinetic model, describe the multi-phase or the multi-component 
dynamics via a diffuse interface approach. There exist therefore a natural regularizing mi- 
croscopic length which is of the order of the interface width, typically a few grid points. Such 
a length scale is also of the order of the effective slip length that must be used whenever a 
quantitative comparison between the hydrodynamical behavior of the LB and the evolution 
of the equivalent Navier-Stokes system is made, as shown for example recently in 15 1. 



IV. COMPARING THE LUBRICATION MODEL AND LATTICE BOLTZMANN 
SIMULATIONS 

In this section we compare the results of the GL model and the LB simulations. Since the 
latter are limited to moderate viscosity ratios i?, the comparison is done for R = 0.03, 0.8 and 
1. We further explore the parameter space in Sec. |Vl using only the lubrication approach. 

The results for this section are computed for 6^ = 7r/2 and a = tt/2. The lattice separation 
in LB is 0.01 (in capillary length units), which will be related to an effective slip length from 
the comparison with the lubrication model. 
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A. Meniscus rise 



We first compare the meniscus rise A for viscosity ratios R = 0.03, 0.8 and 1 in Fig. |31 
Wfien tlie plate is at rest, Ca = 0, tlie meniscus is perfectly horizontal A = due to the 
choice of 6e and a. Let us first consider the case where both liquids have identical viscosity, 
R = 1, but are still immiscible due to the nonzero surface tension. This case is perfectly 
symmetric in the sense that Ca — ?> — Ca gives A —A: there is no difference between 
plunging and withdrawing. This symmetry is indeed observed in Fig. |31 Blue diamonds 
represent LB simulations, while the dash-dotted line is the GL model. We use this symmetric 
case to calibrate the microscopic parameter of the GL model. A nearly perfect fit is achieved 
for slip length = 0.002, which is a reasonable value given that the grid size used in the 
LB simulation is O.OL As Ca increases, the viscous forces increasingly deform the interface, 
leading to a change in A. 

It is interesting to see to what extent the same microscopic parameter is able to describe 
different viscosity ratios. We first mildly decrease the viscosity ratio to R = 0.83 and still 
find a very good agreement between LB an GL (red circles and dotted line respectively). 
With respect to the case R = Iwe see that | A| is slightly smaller at a given value of Ca. This 
means that for the same speed, the meniscus is deformed by a smaller amount due to the 
lower viscosity of the upper phase. When further decreasing the viscosity ratio to R = 0.03 
(green squares, dashed line), some differences between LB and GL becomes apparent (the 
same value for As is maintained). The meniscus in GL is systematically below the value 
obtained in LB. A possible explanation for this difference is the sensitivity of the result on 
the microscopic contact angle imposed as a boundary condition. Still, both models agree 
reasonably well and display very similar trends. In particular, we find that much larger 
values of Ca can be achieved due to the strong reduction of the viscosity in the upper phase. 
This effect is most pronounced for the plunging case, for which the liquid is advancing. This 
is consistent with experimental observations that advancing contact lines can move much 
more rapidly than receding contact lines Q, [s], [sj . The breakdown of the steady solutions, 
which signals the transition to air/liquid entrainment, will be discussed in details in Sec. |Vl 
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FIG. 3: (Color online) Meniscus deformation A as a function of Ca for Be = 7r/2. Ca is nega- 
tive/positive when the plate is moving upward/downward. Symbols: results of lattice Boltzmann 
(LB) simulation. Lattice separation = 0.01. Curves: results of the lubrication- type (GL) model 
with a slip length A^, = 0.002. All lengths are scaled by the capillary length i^. 

B. Shape of the meniscus 

A much more detailed test for the two models is to investigate the detailed structure 
of the interface: How well do the shapes of the menisci compare between GL and LB? In 
Fig. m we show the dynamical meniscus profiles for Ca = 0.019, 0.028 and 0.036, in the 
case of equal viscosities, R = 1. Note that the contact line position is held constant at 
X = 0, so that the bath appears at different heights due to the increase in magnitude of 
A with speed. The agreement of the results of GL model and LB simulation is very good 
in particular for Ca=0.019 and 0.028, even down to the contact line region (Fig. Wp)- For 
larger plate velocities some differences become apparent. These differences may be due to 
different reasons. First, one has to notice that a large Ca is also accompained by a large 
viscous stress contribution, leading to a larger bending of the interface and therefore to the 
possibility to leave the realm of application of the GL model. Second, as said, in the LB 
approach the effective slip length is an output and not an input as for the GL, and it is not 
clear that one would need a finer tuning of it as a function of the capillary number, in order 
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♦ LB, R = 1.0 
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FIG. 4: (Color online) Dynamical meniscus profiles —x vs h for i? = 1 {6e = vr/2 and = 0.002 
for the GL model, lattice separation = 0.01 for the LB simulation). All lengths have been scaled by 
the capillary length. The contact line is at x = so that the bath is at different x for different Ca. 
(a) Solid curves: results of GL model. Dots: results of LB simulation, (b) Zoom on the contact 
line region. 
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to match the GL behavior. 

An even more detailed characterization of the meniscus shape is provided by the local 
angle of the interface Q vs see Fig. O Clearly, both the GL model and the LB sim- 
ulation exhibit the same nontrivial variation of the contact angle. At small scales, the 
angle approaches = 7r/2, while at large scale the meniscus evolves towards the reservoir 
6 = a = 71/2. In between, the angle changes significantly due to the well-known effect of 
"viscous bending" the balance of viscosity and surface tension leads to a curvature of the 
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FIG. 5: (Color online) Local meniscus angle 6 vs h for R = 1. Identical parameters as in Fig. 01 



interface. Very similar variations of the meniscus angle have been obtained experimentally 
40|. 



V. MAXIMUM SPEEDS AND T^NSITION TO AIR ENTRAINMENT 

In this section we discuss the physics of air entrainment in the case of a plunging plate 
(Ca > 0). For realistic situations the viscosity ratio R is typically very small, of order 10~^ 
for water and much smaller for very viscous liquids. This regime can be accessed by the 
GL model only, as LB is restricted to moderate viscosity contrasts. In the first part of this 
section we discuss how the transition to air entrainment is captured in our model in terms 
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Ca =0.0182 




FIG. 6: Meniscus fall A versus Ca {0e = 2.8 radians, R = 0.01, = 10"^). The horizontal dashed 
line indicates the minimum value of A for a static meniscus (= -y/2, with 6e = TT). 



of a bifurcation diagram. Next we study the effect of viscosity ratio on the critical speed. In 
the last part we investigate how the critical speed depends on the microscopic parameters 
such as the slip length and the static contact angle. 



A. Maximum speed for advancing contact lines 



We first consider a case where the equilibrium contact angle is close to vr. For such a 
hydrophobic substrate we expect air entrainment to occur at relatively small Ca jsl, |loi , and 
therefore relatively weak curvature of the interface. This is important, since the assumption 
underlying the GL model is that the interface curvature is weak 28|. We will therefore 
focus on Of. = 2.8 radians and explicitly verify how strongly the interface is curved for our 



numerical solutions. Figure E] shows the drop of the meniscus A as function of Ca for a 
viscosity ratio R = 0.01 (As = 10"^, i.e. of the order of 10 nm). As Ca increases, the contact 
line equilibrates at a lower position resulting in a more negative value of A. However, when 
Ca achieves a certain critical value, stationary solutions cease to exist. This corresponds 
to a maximum plate velocity, or critical capillary number Cac. By analogy to deposition 
of liquid films for plate withdrawal |20l-l22l. l4ll . |42|, this can be identified as the onset of 
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FIG. 7: (Color online) Meniscus profiles for the upper branch (black solid curve) and the lower 
branch (red dashed curve) solutions for Ca = 0.017 {6^ = 2.8 radians, R = 0.01, = 10"^). Here 
we define z = A — x. So the bath level is at z = 0. 



air entrainment: above Cac, unsteady solutions will develop, with a downward motion of 
the contact line js]. As can be seen in Fig. El the critical point arises close to A = —a/2, 
which according to ([T]) corresponds to a meniscus with apparent contact angle tt. This is 
the analogue of the withdrawal case, for which A = +-\/2 and the apparent contact angle 



vanishes at the transition 
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21| . Note that viscous effects push system slightly below 



this maximum extent of deformation for a perfectly static meniscus, with the critical point 
slightly below A = —a/2. 

Interestingly, for a range of speeds Ca < Cac one actually finds more than one solution (cf. 
Fig. [H]). Upon decreasing A, the capillary number first increases and then decreases close 
to the critical point. We refer to the solution branches around Cac as the upper and lower 
branch respectively. Once again, an identical bifurcation structure was previously observed 



for the withdrawing plate case 
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43 



44j . To compare these two types of solutions, we 



show the corresponding meniscus profiles for Ca = 0.017 in Fig. [71 At a large distance from 
the contact line, the solutions are almost identical in shape. Zooming in on the contact line 
region, however, we see the lower branch (red dashed curve) solution displays a "finger" that 
explains the larger magnitude of A with respect to the upper branch. 
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FIG. 8: (Color online) Meniscus fall A versus Ca for different viscosity ratios {9e = 2.8 radians, 
Xs = 10~^). For the case the gas phase has no viscosity, R = 0, steady-state menisci can be 
maintained to arbitrarily large velocity (within our numerical resolution). 
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FIG. 9: (Color online) Critical capillary number Cac versus viscosity ratio R for different slip 
lengths Xs (9^ = 2.8 radians). The dashed line indicates a power law with exponent -1, which is 
valid for large R. 
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B. Effect of viscosity ratio 



A key parameter for the transition to air entrainment is the viscosity ratio R. Figure [S] 
shows the meniscus fall A as function of Ca for different viscosity ratios: i? = 0, 10~^, 
10^'^ and 10~^. At Ca = 0, all cases have the same value of A = —1.15 corresponding to a 
static bath with contact angle 9^, = 2.8. Now we consider the result for R = 0, for which 
there are no viscous effects in the air {r]g = 0). We observe that A decreases with Ca, but 
without any bifurcation. It appears that steady meniscus solutions can be sustained up to 
arbitrarily large plate velocities. In fact, the curve is consistent with the scaling A ~ — y^Ca 
at relatively large Ca, corresponding to a simple balance between gravity and viscosity. For 
R ^ 0, however, the situation becomes fundamentally different. While the curves follow the 
same trend as for i? = at small Ca, a deviation appears at larger speeds that ultimately 
leads to a critical point. Each nonzero viscosity ratio has a well-defined critical speed, with 
Cac increasing when the viscosity ratio R is reduced. 

These observations can be interpreted as follows. As long as the viscosity of the air has a 
negligible effect on the flow, the curves are indistinguishable from the case R = 0. Deviations 
of the -R = curve signal that the air flow starts to influence the shape of the meniscus. 
Physically, this arises because the interface slope approaches vr, leaving only a narrow wedge 
angle for the air flow. Figure [2] illustrates that the recirculation in the air induces significant 
velocity gradients: despite the small air viscosity, the stresses in the small wedge of air 
become comparable to those in the liquid. Mathematically, this can be derived from the 
function f{6, R) as defined in f|T6|) . For small R and 6 close to vr, we can approximate: 

f{9, R) ^ f{9, 0) - 4i? ^ _2{n_el_ _ ^^^^ 

ovr 

as long as vr — ^ ^ 27rR. For 6 very close to vr, f{6, R) has a different asymptotic form, see 
the appendix for details. The first term in ( 13T]) represents the (relative) viscous contribution 
inside the liquid, which vanishes in the limit ^ — t- tt. The second term represents the viscous 
contribution in the air, which will be significant once (tt — ^) ~ R^^^. Noting that the contact 
line is receding from the point of view of the air phase, one understands that a critical speed 
appears when the effect of the air becomes important. 

The dependence of the critical speed Cac on the viscosity ratio R is shown in Fig. [9] 
(for various slip lengths). First we consider the limit R ^ 1, for which the upper fluid 
is actually much more viscous than the bottom fluid. This is the usual case of a receding 
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contact line that is completely dominated by the upper (receding) phase. In this limit we 
expect the critical speed to scale with the viscosity of the upper phase, denoted 77^, such 
that Uc ~ l/Vg- Since we have based the capillary number on the viscosity rjg, we obtain 
Cac = UcTji/^ ~ . This is indeed observed in Fig. |9] at i? ^ 1. However, our main 
interest here lies in the opposite limit, i.e. i? ^ 1, as for air entrainment. As already 
mentioned, the critical speed seems to increase indefinitely by reducing the viscosity ratio. 
This suggests that for the limiting case of i? = 0, steady menisci can be sustained at 
arbitrarily large plunging velocities. Our numerical resolution does not allow for a perfect 
determination of the asymptotics for i? ^ 1. Enforcing a power law fit, Cac ~ -R^, in the 
range R = 10~^ — 10^^, one obtains /3 = —0.67. This (effective) exponent suggests that 
both phases play an important role in determining the critical speed. Namely, the exponent 
would be /3 = — 1 if only r]g were important, while /3 = corresponds to the case where rj^ 
is the only relevant viscosity. 

Finally, we briefly verify the assumption of small curvature, necessary for the strict va- 
lidity of the model. In Fig. [10] we plot the dimensionless curvature, h\d6/ds\, as function of 
h in the vicinity of the critical point {R = 10"^, 10"^ and 10"^, = 10"^). At small scales, 
h\d6/ds\ ^ 1 for all Ca, consistent with the assumption of small curvature. However, the 
curvature increases significantly when approaching the bath due to the bending of interface 
from a large contact angle to tt/2. The magnitude is acceptable in this regime, in particular 
since viscous effects becomes less important at large scales. Inclining the plate angle to 
values close to n would further reduce this bending effect, and extend the range of vahdity 
of the GL model. 

C. Dependence of the critical speed on microscopic parameters 

Apart from the viscosity ratio, the GL model contains two parameters: the slip length 
and the microscopic (equilibrium) contact angle Of.. Here we discuss the dependence of Cac 
on these parameters. The slip length was varied already in Fig. [HI with values A^ = 10~^, 
10~^ and 10~^. As expected for wetting problems, we see a weak increase of Cac with 
Xg. A larger A^ reduces the range over which viscous dissipation is effective. This leads 
to a (logarithmic) reduction of the viscous dissipation, while the capillary driving remains 
unaltered. Hence, larger velocities can be achieved before air entrainment occurs. 
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FIG. 10: Scaled curvature \hd9/ds\ vs h for Ca very close to Cac (A^ = 10~^). 

The dependence of the critical speed on Of, is investigated in Fig. [TT] The figure reveals 
that there is no obvious universal scaling behavior for Cac down to viscosity ratios as small as 
R = 10~^. Enforcing a power-law fit, different 6^ would give rise to different exponents. We 
do clearly see that Cac decreases with 6e, which is further emphasized in Fig. [T2j Consistent 
with jo], Ilo|, the critical speed vanishes in the hydrophobic limit where Of, — ?■ vr. We note 
that for contact angles that are not close to vr, the shape of the meniscus displays significant 
curvatures. In this sense, we expect that our results are not fully quantitative solutions of 
the Stokes fiow problem. 

VI. DISCUSSION AND CONCLUSIONS 

In this paper we present, compare and employ two distinct models to study the meniscus 
deformation and the onset of air entrainment in a dip-coating geometry. The first model is a 
generalization of the lubrication theory to a two-phase fiow situation, in which a slip length 
is introduced to resolve the viscous singularity. The second model is a numerical one based 
on the discretization of the Boltzmann equation, namely the Lattice Boltzmann method for 
multiphase/multicomponent fiuids. In this model the viscous singularity is removed by the 
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FIG. 11: (Color online) Critical speed Cac as a function of R for different static contact angle 9e 
(As = 0.001). Symbols are results of GL model. Dashed curves are predictions of Cox's model (top 
one for 9(. = 0.87 radian, bottom one for 6(, = 3.0 radian), from Eq. (8.3) of for which we take 
the ratio between the microscopic length scale and macroscopic length scale to be 10~^ (the same 
value as our slip length scaled by the capillary length A^). The dashed-dotted line indicates slope 
of -1. 



lattice discretization, which effectively introduces an effective slip length to the system [45 1. 
The effective slip length also depends on the choice of multiphase/multicomponent model. 

The results of GL and LB have a good agreement, in particular when Ca is relatively 
small. When exploring larger values of Ca, the two models start to differ as shown in Fig. 
ini which can be attributed to different physics at microscopic and hydrodynamical scales 
(e.g. the non-zero interface thickness in the LB or strong viscous bending of interface breaks 
the hypothesis of small interface curvatures made in the GL approach). Yet, qualitative 
features, such as the bending of the meniscus and the dependence on viscosity ratio, are 
consistent for the two models. The transition to air entrainment for 6^ close to vr involves 
relatively weak curvatures and is thus captured by the GL model. For the LB simulations 
the main challenge is given by the large viscosity contrasts, which is still not fully achievable 
for the multiphase/multicomponent LB model used here. 
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FIG. 12: Critical speed Cac as a function of static contact angle 9e for R = 0.001 and = 0.001. 

In the second part of this paper, the critical speed of air entrainment is investigated by the 
GL model. We have found a strong dependence of critical speeds on the air viscosity, which 
is consistent with the experiments performed by Marchand et al. js]. Remarkably, both our 
theoretical results and the experimental results from js] differ from Cox's model in which 



Cac is predicted to depend only logarithmically on air viscosity 12[. For comparison, we 
have added in Fig. [11] predictions of Cox's model, from Eq. (8.3) of 12|, represented by the 
dashed curves (top curve for = 0.87 radian, bottom one for O^, = 3.0 radian). For large R, 
we find that both models predict Cac scales as R with scaling —1. This regime corresponds 
to the usual dewetting case for which the critical speed only depends on the viscosity of the 



more viscous fluid [17 
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22| . Note that in Cox's model, there is an undetermined factor 



e which is defined as the ratio of the microscopic length scale to the macroscopic length 
scale. Here we take e to be the same value as As, which is 10~^. Interestingly, both the 
models predict exactly the same values of Cac for 6^ = 0.87 radian when R is large, which 
we consider as an coincidence since e is an adjustable parameter. For 6e = 3.0 radian. Cox's 
and our results differ by a factor. More interesting things occur at small R, it is clearly 
shown in Fig. [TT] that for Cox's model (red top curve ) Cac increases extremely slowly 
(logarithmically) as R is decreased. By contrast, our model predicts a moderate increase of 
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Cac- Interestingly, such a weak logarithmic relation has been observed in the case of liquid 
impacting on liquid [isl [l^, for which there is no moving contact line. Both our theoretical 
results and experimental results from 8[ therefore suggest that the mechanism leading to air 
entrainment can be fundamentally different depending on whether a contact line is present 
or not. 

The GL model is directly compared with experiments in js], and shows that the model is 
able to qualitatively capture the dependence of Cac on the viscosity ratio rjg/rji. Quantita- 
tively, however, the agreement is not satisfactory (see Fig. 3 of js]). We believe this is due 
to the relatively large meniscus curvatures encountered in the experiments (static contact 
angle of the substrate ~ 50°), pushing the problem beyond the assumptions of the model. 
It would be interesting to explore other methods to achieve a more quantitative description 
of air entrainment by advancing contact line, in particular for large values of Ca. From an 
experimental perspective, more insight could be obtained by varying the gas viscosity or 
by replacing the air by a liquid of low viscosity. It would also be interesting to perform 
experiments with a substrate of large static contact angles. 
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Appendix A: Appendix 

We discuss the behavior of the function f{6, R) when 6 is close to vr and R close to zero. 
As presented in Section UlI A \ \ f{9,R) is defined as: 



f{0,R) 



2sm^e[R'Me) + 2Rfsie) + h{7r- 

3[Rf,(e)M7r-e)-f,{n-e)f2 



25 



fi{9) = e^-siii^e 

f2{9) = ^ -sin^cos^ 
fs{e) = {9(71 -9) + sin' 9). 

(Al) 

First, we expand the terms in both the numerator and the denominator in series of (vr — ^) 
and keep the leading order terms only, we end up with 

The asymptotic behavior of /(6', R) depends on the relative magnitude between (tt — 9) 
and R. For R <^ -k — 9, f{9,R) can be approximated as: 

fi9, R) f{9, 0) - 4i? -'^(^ -(^)' _ (A3) 

ovr 

The contribution of air viscosity, represented by —4R, will become significant once {71 — 9) ~ 
R'/\ 

When 9 is very close to vr such that tt — 9 R, f {9 , R) goes to a different asymptotic 
form, i.e. 

f{9,R)c^-R. (A4) 

If we substitute this asymptotic form of f{9,R) into the generalized lubrication equation 
( ]22|) . we will see the liquid viscosity will be canceled out in the multiple CaR so that the 
asymptotic equation does not depend on liquid viscosity. This means in this asymptotic 
limit, air viscosity completely dominates the flow. 
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